function ImageBinary=plotVisibility_both(IC,T,u,FN_i,FN_r,CameraAngleCapability,startTime,DU,VU)

options=odeset('RelTol',1e-13,'AbsTol',1e-13);
%% Propagate state in rotating and calculate state in inertial
to = [0:0.01:T];
Xo = IC;
[t,S] = ode113(@(t,S)CR3BP_n(t,S,u),to,Xo,options);
S = S';
S_i = zeros(size(S));
for ii = 1:length(S)
    S_i(:,ii) = C2I_primary(S(:,ii),u,DU,VU,t(ii)+startTime);
end
%% Inertial

% plot3(S(1,:),S(2,:),S(3,:))
[row,col,~] = size(FN_i);
visibility_i=zeros(row,col);
parfor ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = 1:length(t)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if angle_temp < CameraAngleCapability
                visibility_i(ii,jj) = 1;
            end
        end
    end
end


ImageBinary = mat2gray(visibility_i);

Window = figure('OuterPosition',[100,100,1500,800]); hold on;
mesh(ImageBinary*100,'EdgeColor','none','FaceColor','flat');
xlabel('Longitude')
ylabel('Latitude')
title('Sun Visibility (inertial frame)')
set(gca,'fontsize',16)
xticks(linspace(0,row,5));
xticklabels({'180W','90W','0','90E','180E'});
yticks(linspace(0,col,5))
yticklabels({'90S','45S','0','45N','90N'});
xlim([1,row]);
ylim([1,col])

%% Inertial

% plot3(S(1,:),S(2,:),S(3,:))
% [row,col,~] = size(FN_r);
% visibility_r=zeros(row,col);
% parfor ii = 1:row
%     for jj = 1:col
%         normvec = [FN_r(ii,jj,1);FN_r(ii,jj,2);FN_r(ii,jj,3)];
%         for k = 1:length(t)
%             posvec = S(1:3,k);
%             angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
%             if angle_temp < CameraAngleCapability
%                 visibility_r(ii,jj) = 1;
%             end
%         end
%     end
% end
% 
% 
% ImageBinary = mat2gray(visibility_r);
% 
% Window = figure('OuterPosition',[100,100,1500,800]); hold on;
% mesh(ImageBinary*100,'EdgeColor','none','FaceColor','flat');
% xlabel('Longitude')
% ylabel('Latitude')
% title('Sun Visibility (rotating frame)')
% set(gca,'fontsize',16)
% xticks(linspace(0,row,5));
% xticklabels({'180W','90W','0','90E','180E'});
% yticks(linspace(0,col,5))
% yticklabels({'90S','45S','0','45N','90N'});
% xlim([1,row]);
% ylim([1,col])
% a=colorbar;
% ylabel(a,'Access %','FontSize',16,'Rotation',270);
% a.Label.Position(1) = 4;
end